### Supplementary Materials for the Article "Rethinking Public Funding of Parties and Corruption: Confronting Theoretical Complexity and Challenging Measurement"

#### The code from this script reproduces all tables and figures from the online supplementary appendix for the article "Rethinking Public Funding of Parties and Corruption: Confronting Theoretical Complexity and Challenging Measurement", with the exception of Table A1 that presents data sources


### Set your working directory
setwd("your_working_directory")

### Install the required packages if you do not have them installed 
#install.packages(tidyverse)
#install.packages(haven)
#install.packages(fixest)
#install.packages(geepack)
#install.packages(panelr)
#install.packages(broom)
#install.packages(modelsummary)
#install.packages(bestNormalize)
#install.packages(ggeffects)
#install.packages(cowplot)
#install.packages(scales)
#install.packages(corrplot)
#install.packages(ggrepel)
#install.packages(kableExtra)

### Download libraries 

library(tidyverse)
library(haven)
library(fixest)
library(geepack)
library(panelr)
library(broom)
library(modelsummary)
library(bestNormalize)
library(ggeffects)
library(cowplot)
library(scales)
library(corrplot)
library(ggrepel)
library(kableExtra)

### Import and manipulate data 

impact <- read_dta("impact.dta") %>%
  mutate(year_factor = factor(year), 
         country_factor = factor(country))
         
################################################
### Table of summary statistics

### Variable selection and renaming
data_summary <- impact %>%
  select(country, capture_impact_1, capture_impact_2, capture_impact_3, capture_impact_1_norm, capture_impact_2_norm, capture_impact_3_norm, dpfv_total, dpfv_nomin_infl, dpfv_dprc_infl, private_agreg, democracy_scaled, gdppc_ppp_2017, elect_syst, govern_type2, gini_disp, gov_size_imf, ethn_fr) %>%
  mutate(gdppc_ppp_2017 = gdppc_ppp_2017/1000) %>%
  rename("Country" = "country", 
         "Corruption 1" =  "capture_impact_1", 
         "Corruption 2" =  "capture_impact_2", 
         "Corruption 3" =  "capture_impact_3",
         "Corruption 1 Yeo-Johnson transformed" =  "capture_impact_1_norm", 
         "Corruption 2 Yeo-Johnson transformed" =  "capture_impact_2_norm", 
         "Corruption 3 Yeo-Johnson transformed" =  "capture_impact_3_norm",
         "DPF level nominal value" = "dpfv_total", 
         "DPF level inflation adjusted" = "dpfv_nomin_infl", 
         "DPF level inflation adjusted & discounted" = "dpfv_dprc_infl",
         "Donation limits restrictions" = "private_agreg", 
         "Democracy" = "democracy_scaled",
         "GDP per capita (.000)" =  "gdppc_ppp_2017", 
         "Electoral system" = "elect_syst", 
         "Parliamentary" = "govern_type2",
         "Inequality" = "gini_disp",
         "Government size % GDP" = "gov_size_imf",
         "Ethnic fractionalization"  = "ethn_fr")

### Saummary table for appendix A
sum_data <- datasummary(`Corruption 1` + `Corruption 2` + `Corruption 3` + `Corruption 1 Yeo-Johnson transformed` + `Corruption 2 Yeo-Johnson transformed` + `Corruption 3 Yeo-Johnson transformed` + `DPF level nominal value` + `DPF level inflation adjusted` + `DPF level inflation adjusted & discounted` + `Donation limits restrictions` + Democracy + `GDP per capita (.000)` + `Electoral system` + Parliamentary + Inequality + `Ethnic fractionalization` + `Government size % GDP` ~ Mean + Median + Min + Max + SD + P25 + P75, data = data_summary, output = 'data.frame') 
#########################################
### Table A2 online appendix
kbl(sum_data, booktabs = TRUE, longtable = T, caption = "Table of summary statistics") %>%
  kable_styling(font_size = 12, full_width = F, latex_options = c("striped", "hold_position")) %>%
  column_spec(column = 1, width = "3.5cm") %>%
  column_spec(column = 2, width = "1.5cm") %>% 
  column_spec(column = 3, width = "1.5cm") %>% 
  column_spec(column = 4, width = "1.5cm") %>% 
  column_spec(column = 5, width = "1.5cm") %>%
  column_spec(column = 6, width = "1.5cm") %>% 
  column_spec(column = 7, width = "1.5cm") %>% 
  column_spec(column = 8, width = "1.5cm") %>%
  footnote(general = "Corruption 1 - percentage of informal payments to political parties and parliamentarians; Corruption 2 - percentage of informal payments to political parties, parliamentarians and central government officials; Corruption 3 - percentage of informal payments to political parties, parliametarians, central and local government officials. \\\\textit{Corruption 1 Yeo-Johnson transformed} and \\\\textit{DPF level inflation-adjusted and inflation-adjusted and discounted} are used as dependent and and independent variables in the main analysis respectively.", escape = FALSE, fixed_small_size = TRUE, footnote_as_chunk = TRUE, threeparttable = TRUE)

###############################################
### Correlation matrix for the appendix A

summary_matrix <- cor(data_summary[, c(2, 5, 8:18)]) 

### Plot and export the matrix, Figure A1 appendix

corrplot(summary_matrix,
         method = "square",
         order = "hclust",
         tl.col = "grey20", 
         tl.srt = 90, 
         number.cex = 0.8,
         addCoef.col = "grey40",
         diag=FALSE,
         type = "upper")

##############################################
### Figure A2 appendix, Relationship between alternative operationalisations of public funding

p1 <- ggplot(impact, aes(x = dpfv_total, y = dpfv_nomin_infl)) +
  geom_point(shape = 1, size = 2, col = "lightskyblue4") +
  geom_smooth(method = "lm", color = "lightsteelblue4") +
  scale_y_continuous(breaks = seq(0,8,2)) +
  theme_light(base_family = "serif") +
  labs(x = "DPF nominal value", y = "DPF inflation adjusted")

p2 <- ggplot(impact, aes(x = dpfv_total, y = dpfv_dprc_infl)) +
  geom_point(shape = 1, size = 2, col = "lightskyblue4") +
  geom_smooth(method = "lm", color = "lightsteelblue4") +
  scale_y_continuous(breaks = seq(0,8,2)) +
  theme_light(base_family = "serif") +
  labs(x = "DPF nominal value", y = "DPF inflation adjusted & discounted")

p3 <- ggplot(impact, aes(x = dpfv_nomin_infl, y = dpfv_dprc_infl)) +
  geom_point(shape = 1, size = 2, col = "lightskyblue4") +
  geom_smooth(method = "lm", color = "lightsteelblue4") +
  scale_y_continuous(breaks = seq(0,8,2)) +
  theme_light(base_family = "serif") +
  labs(x = "DPF inflation adjusted", y = "DPF inflation adjusted & discounted")

cowplot::plot_grid(p1, p2, p3, ncol = 3, align = "h")

#########################################
### Figure A3 appendix, Relationship between public funding per registered voter and vote

p4 <- ggplot(impact, aes(x = dpfrv_total + 0.01, y = dpfv_total + 0.01)) +
  geom_point(size = 2, shape = 1, color = "lightskyblue4")+
  geom_smooth(method = "lm", color = "lightsteelblue4") +
  scale_x_log10() +
  scale_y_log10() +
  theme_light(base_family = "serif") +
  labs(x = "DPF per registered voter: \nnominal value (log10 axis)",
         y = "DPF per vote: \nnominal value (log10 axis)")

p5 <- ggplot(impact, aes(x = dpfrv_total_infl + 0.01, y = dpfv_nomin_infl + 0.01)) +
  geom_point(size = 2, shape = 1, color = "lightskyblue4")+
  geom_smooth(method = "lm", color = "lightsteelblue4") +
  scale_x_log10() +
  scale_y_log10() +
  theme_light(base_family = "serif") +
  labs(x = "DPF per registered voter: \ninflation-adjusted (log10 axis)",
         y = "DPF per vote: \ninflation-adjusted (log10 axis)")

cowplot::plot_grid(p4, p5, ncol = 2, align = "h")

###########################################
### Figure A4 appendix, Relationship between alternative operationalisation of corruption
p6 <- ggplot(impact, aes(x = capture_impact_1, y = capture_impact_2)) +
  geom_point(shape = 1, size = 2, col = "lightskyblue4") +
  geom_smooth(method = "lm", color = "lightsteelblue4") +
  theme_light(base_family = "serif") +
  labs(x = "Party & legislative corruption", y = "Party, legislative & central \ngovernment corruption")

p7 <- ggplot(impact, aes(x = capture_impact_1, y = capture_impact_3)) +
  geom_point(shape = 1, size = 2, col = "lightskyblue4") +
  geom_smooth(method = "lm", color = "lightsteelblue4") +
  theme_light(base_family = "serif") +
  labs(x = "Party & legislative corruption", y = "Party, legislative, central & \nlocal government corruption")

cowplot::plot_grid(p6, p7, ncol = 2, align = "h")

#######################################
## Normalization of the dependent variable

set.seed(546)

### Alternative normalizations of the dependent variable 

corruption_normalized <- bestNormalize::bestNormalize(impact$capture_impact_1) # Normaization 

corruption_normalized # Output with the selection of the best transformation

# Extract the best transformation
best_df <- map(corruption_normalized$chosen_transform, ~unlist(.)) %>%
  bind_cols() %>%
  select(2:1) %>%
  rename("original_values" = "x", "yeo_johnson" = "x.t")

# Extract other transformations
best_df2 <- map(corruption_normalized$other_transform, ~unlist(.[[1]])) %>%
  bind_cols()

### Bind all transformations together and rename
all_transformations <- cbind(best_df, best_df2) %>%
  pivot_longer(cols = c(1:10), names_to = "Transformations", values_to = "Score") %>%
  mutate(transformations_2 = 
           case_when(Transformations == "original_values" ~ "Original values",
                     Transformations == "yeo_johnson" ~ "Yeo-Johnson", 
                     Transformations == "arcsinh_x" ~ "Arcsinh",
                     Transformations == "boxcox" ~ "Box-Cox",
                     Transformations == "center_scale" ~ "Center-Scale",
                     Transformations == "double_reverse_log" ~ "Double Reversed Log",
                     Transformations == "exp_x" ~ "Exponential",
                     Transformations == "log_x" ~ "Log",
                     Transformations == "orderNorm" ~ "Ordered-quantile",
                     Transformations == "sqrt_x" ~ "Sqrt"))

### Plot transformations, Figure A5 appendix

all_transformations %>%
  filter(!Transformations %in% c("exp_x", "double_reverse_log")) %>%
  ggplot(aes(x = Score, color = transformations_2)) +
  geom_density() + 
  facet_wrap(~ transformations_2, scales = "free", ncol = 2) +
  theme_minimal(base_family = "serif") +
  theme(legend.position = "null",
        panel.grid.minor = element_blank()) +
  labs(x = "", y = "")

############################################

### Fit year fixed effects with public funding inflation adjusted and discounted 
### Different combinations of controls
### Standard errors clustered by country

### Declare panel data 
impact_panel <- panel(impact, ~country + year)

### Fit Year FE models 
models_yfe <- list(
  m1_yfe <- feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | year, data = impact_panel, cluster = ~country), ## Model 1
  update(m1_yfe, . ~ . + elect_syst), # Model 2
  update(m1_yfe, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_yfe, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_yfe, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_yfe, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

## Clean term labels for the regression table
model_labels_yfe <- c("dpfv_dprc_infl" = "DPF level", 
                  "private_agreg" = "Donation limits", 
                  "democracy_scaled" =  "Democracy", 
                  "log(gdppc_ppp_2017)" = "GDPpc(log)",
                  "elect_syst" = "Electoral System", 
                  "govern_type2" = "Parliamentary",
                  "ethn_fr" = "Ethnic fraction.", 
                  "gini_disp" = "Inequality", 
                  "gov_size_imf" = "Government size % GDP")

# Export regression table, Table B1 appendix
modelsummary(models_yfe,
             stars = TRUE,
             coef_map = model_labels_yfe,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Year fixed-effects models (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

#########################################
### Fit country fixed effects with public funding inflation adjusted and discounted
### Different combinations of controls
### Standard errors clustered by country

### Fit country FE
models_cfe <- list(
  m1_cfe <- feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country, data = impact_panel, cluster = ~country), # Model 1
  update(m1_cfe, . ~ . + elect_syst), # Model 2
  update(m1_cfe, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_cfe, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_cfe, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_cfe, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

# Export regression table, Table B2 appendix
modelsummary(models_cfe,
             stars = TRUE,
             coef_map = model_labels_yfe,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Country fixed-effects models (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")


############################################
### Fit two-way fixed effects with public funding inflation adjusted and discounted
### Different combinations of controls
### Standard errors clustered by country
### Fit TWFE
models_tfe <- list(
  m1_tfe <- feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country + year, data = impact_panel, cluster = ~ country), # Model 1
  update(m1_tfe, . ~ . + elect_syst), # Model 2, 
  update(m1_tfe, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

# Export regression table, Table B3 appendix
modelsummary(models_tfe,
             stars = TRUE,
             coef_map = model_labels_yfe,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Two-way fixed-effects models (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

#########################################
### Fit Within-Between Random-Effects (WBRE) models  with public funding inflation adjusted and discounted

#### Declaring panel data
panel_wbre <- panel_data(impact, id = country, wave = year) 
### Fit WBRE
wbre_models <- list(
  m1_wbgee <- wbgee(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017), data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1"), # model 1
  update(m1_wbgee, .~. + elect_syst), # model 2
  update(m1_wbgee, .~. + elect_syst + govern_type2), # model 3
  update(m1_wbgee, .~. + elect_syst + govern_type2 + gini_disp), # model 4
  update(m1_wbgee, .~. + elect_syst + govern_type2 + ethn_fr), # model 5
  update(m1_wbgee, .~. + elect_syst + govern_type2 + gov_size_imf)) # model 6

# Clean term labels labels
coeff_labels_wbre <- 
  c("dpfv_dprc_infl" = "DPF level (within)", 
    "private_agreg" = "Donation limits (within)", 
    "democracy_scaled" =  "Democracy (within)", 
    "log(gdppc_ppp_2017)" = "GDPpc(log) (within)", 
    "elect_syst" = "Electoral System (within)", 
    "govern_type2" = "Parliamentary (within)",
    "ethn_fr" = "Ethnic fractionalization (within)", 
    "gini_disp" = "Inequality (within)", 
    "gov_size_imf" = "Government size % GDP (within)", 
    "imean(dpfv_dprc_infl)" = "DPF level (between)", 
    "imean(private_agreg)" = "Donation limits (between)", 
    "imean(democracy_scaled)" =  "Democracy (between)",
    "imean(log(gdppc_ppp_2017))" = "GDPpc(log) (between)", 
    "imean(elect_syst)" = "Electoral System (between)", 
    "imean(govern_type2)" = "Parliamentary (between)",
    "imean(ethn_fr)" = "Ethnic fractionalization (between)", 
    "imean(gini_disp)" = "Inequality (between)", 
    "imean(gov_size_imf)" = "Government size % GDP (between)")

### Table B4 appendix
modelsummary(wbre_models, 
             stars = TRUE, 
             coef_map = coeff_labels_wbre,
             caption = "The effect of DPF on party corruption: Within-between random-effects models (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             longtable = TRUE,
             output = "kableExtra") %>%
  kable_styling(font_size = 10, full_width = F, latex_options = "hold_position") %>%
  column_spec(column = 1, width = "4cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

##########################################

### Fit year fixed effects with public funding inflation adjusted and discounted
### Different combinations of controls
### Fit Year FE
models_yfe2 <- list(
  m1_yfe2 <- feols(capture_impact_1_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | year, data = impact_panel, cluster = ~country), ## Model 1
  update(m1_yfe2, . ~ . + elect_syst), # Model 2
  update(m1_yfe2, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_yfe2, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_yfe2, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_yfe2, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

## Model labels for regression tables
model_labels_yfe2 <- c("dpfv_nomin_infl" = "DPF level", 
                  "private_agreg" = "Donation limits", 
                  "democracy_scaled" =  "Democracy", 
                  "log(gdppc_ppp_2017)" = "GDPpc(log)",
                  "elect_syst" = "Electoral System", 
                  "govern_type2" = "Parliamentary",
                  "ethn_fr" = "Ethnic fraction.", 
                  "gini_disp" = "Inequality", 
                  "gov_size_imf" = "Government size % GDP")

# Export regression table, Table B5 appendix
modelsummary(models_yfe2,
             stars = TRUE,
             coef_map = model_labels_yfe2,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Year fixed-effects models (DPF inflation adjusted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

########################################
### Fit year fixed effects with public funding inflation adjusted and discounted
### Different combinations of controls
### Standard errors clustered by country
### Fit country FE
models_cfe2 <- list(
  m1_cfe2 <- feols(capture_impact_1_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country, data = impact_panel, cluster = ~country), # Model 1
  update(m1_cfe2, . ~ . + elect_syst), # Model 2
  update(m1_cfe2, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_cfe2, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_cfe2, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_cfe2, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

# Export regression table, Table B6 appendix
modelsummary(models_cfe2,
             stars = TRUE,
             coef_map = model_labels_yfe2,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Country fixed-effects models (DPF inflation adjusted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

########################################
### Fit two-way fixed effects with public funding inflation adjusted and discounted
### Different combinations of controls
### Standard errors clustered by country
### Fit TWFE
models_tfe2 <- list(
  m1_tfe2 <- feols(capture_impact_1_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country + year, data = impact_panel, cluster = ~ country), # Model 1
  update(m1_tfe2, . ~ . + elect_syst), # Model 2, 
  update(m1_tfe2, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_tfe2, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_tfe2, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_tfe2, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

# Export regression table, Table B7 appendix
modelsummary(models_tfe2,
             stars = TRUE,
             coef_map = model_labels_yfe2,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Two-way fixed-effects models (DPF inflation adjusted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

#######################################
### Fit two-way fixed effects with public funding inflation adjusted 
### Different combinations of controls
### Fit WBRE
wbre_models2 <- list(
  m1_wbgee2 <- wbgee(capture_impact_1_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017), data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1"), # model 1
  update(m1_wbgee2, .~. + elect_syst), # model 2
  update(m1_wbgee2, .~. + elect_syst + govern_type2), # model 3
  update(m1_wbgee2, .~. + elect_syst + govern_type2 + gini_disp), # model 4
  update(m1_wbgee2, .~. + elect_syst + govern_type2 + ethn_fr), # model 5
  update(m1_wbgee2, .~. + elect_syst + govern_type2 + gov_size_imf)) # model 6

coeff_labels_wbre2 <- 
  c("dpfv_nomin_infl" = "DPF level (within)", 
    "private_agreg" = "Donation limits (within)", 
    "democracy_scaled" =  "Democracy (within)", 
    "log(gdppc_ppp_2017)" = "GDPpc(log) (within)", 
    "elect_syst" = "Electoral System (within)", 
    "govern_type2" = "Parliamentary (within)",
    "ethn_fr" = "Ethnic fractionalization (within)", 
    "gini_disp" = "Inequality (within)", 
    "gov_size_imf" = "Government size % GDP (within)", 
    "imean(dpfv_nomin_infl)" = "DPF level (between)", 
    "imean(private_agreg)" = "Donation limits (between)", 
    "imean(democracy_scaled)" =  "Democracy (between)",
    "imean(log(gdppc_ppp_2017))" = "GDPpc(log) (between)", 
    "imean(elect_syst)" = "Electoral System (between)", 
    "imean(govern_type2)" = "Parliamentary (between)",
    "imean(ethn_fr)" = "Ethnic fractionalization (between)", 
    "imean(gini_disp)" = "Inequality (between)", 
    "imean(gov_size_imf)" = "Government size % GDP (between)")

### Table B8 appendix
modelsummary(wbre_models2, 
             stars = TRUE, 
             coef_map = coeff_labels_wbre2,
             caption = "The effect of DPF on party corruption: Within-between random-effects models (DPF inflation adjusted)",
             align = "lrrrrrr",
             longtable = TRUE,
             output = "kableExtra") %>%
  kable_styling(font_size = 10, full_width = F, latex_options = "hold_position") %>%
  column_spec(column = 1, width = "4cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")


########################################
### Fit two-way fixed effects with elections dummies, elections frequency and both versions of public funding

models_twfe_dummy <- list(
  
  m1_twfe_d <- feols(capture_impact_1_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf + elections_dummy | country + year, data = impact_panel, cluster = ~ country), # Model 1

    m2_twfe_d <- feols(capture_impact_1_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf + elections_frequency | country + year, data = impact_panel, cluster = ~ country), # Model 2
  
  m3_twfe_d <- feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf + elections_dummy | country + year, data = impact_panel, cluster = ~ country), # Model 3
  
  m4_twfe_d <- feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf + elections_frequency | country + year, data = impact_panel, cluster = ~ country)) # Model 4

# Clean labels 
model_labels_el_dummy <- 
  c("dpfv_nomin_infl" = "DPF level-inflation adjusted",
    "dpfv_dprc_infl" = "DPF level-inflation adjusted & discounted",
    "private_agreg" = "Donation limits",
    "democracy_scaled" =  "Democracy",
    "log(gdppc_ppp_2017)" = "GDPpc(log)",
    "elect_syst" = "Electoral System", 
    "govern_type2" = "Parliamentary",
    "ethn_fr" = "Ethnic fraction.",
    "gini_disp" = "Inequality",
    "gov_size_imf" = "Government size % GDP",
    "elections_dummy" = "Elections dummies",
    "elections_frequency" = "Elections frequency")

## Table B9 appendix
modelsummary(models_twfe_dummy,
             stars = TRUE,
             coef_map = model_labels_el_dummy,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Two-way fixed-effects models with election dummies and election frequency",
             align = "lrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "5.5cm") %>%
  column_spec(column = 2, width = "2cm") %>% 
  column_spec(column = 3, width = "2cm") %>% 
  column_spec(column = 4, width = "2cm") %>% 
  column_spec(column = 5, width = "2cm") 

###########################################
### Fit two-way fixed effects with public funding inflation adjusted and discounted:
### Removing authoritarian regimes
impact_panel_no_authorit <- impact_panel %>%
  filter(!country %in% c("AZE", "BLR", "KAZ", "TJK", "UZB", "KGZ", "RUS"))

### Fit TWFE 
models_tfe_no_authorit <- list(
  m1_tfe <- feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country + year, data = impact_panel_no_authorit, cluster = ~ country), # Model 1
  update(m1_tfe, . ~ . + elect_syst), # Model 2, 
  update(m1_tfe, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

# Export regression table, Table B10 appendix
modelsummary(models_tfe_no_authorit,
             stars = TRUE,
             coef_map = model_labels_yfe,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Two-way fixed-effects models without authoritarian regimes (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country. Authoritarian regimes removed: Azerbaijan, Belarus, Kazakhstan, Kyrgyzstan, Russia, Tajikistan, Uzbekistan", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

########################################
### Fit two-way fixed effects with public funding inflation adjusted
### Authoritarian regimes removed

### Fit TWFE
models_tfe_no_authorit2 <- list(
  m1_tfe <- feols(capture_impact_1_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country + year, data = impact_panel_no_authorit, cluster = ~ country), # Model 1
  update(m1_tfe, . ~ . + elect_syst), # Model 2, 
  update(m1_tfe, . ~ . + elect_syst + govern_type2), # Model 3
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + gini_disp), # Model 4
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + ethn_fr), # Model 5
  update(m1_tfe, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

# Export regression table, Table B11 appendix
modelsummary(models_tfe_no_authorit2,
             stars = TRUE,
             coef_map = model_labels_yfe2,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: Two-way fixed-effects models without authoritarian regimes (DPF inflation adjusted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country. Authoritarian regimes removed: Azerbaijan, Belarus, Kazakhstan, Kyrgyzstan, Russia, Tajikistan, Uzbekistan", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

###########################################
### Fit WBRE models without authoritarian regimes, with public funding inflation adjusted and discounted 

### Declaring panel data
panel_wbre_no_authorit1 <- panel_data(impact, id = country, wave = year) %>% 
  filter(!country %in% c("AZE", "BLR", "KAZ", "TJK", "UZB", "KGZ", "RUS"))
  
### Fit WBRE
wbre_models_no_authorit1 <- list(
  m1_wbgee <- wbgee(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017), data = panel_wbre_no_authorit1, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1"), # model 1
  update(m1_wbgee, .~. + elect_syst), # model 2
  update(m1_wbgee, .~. + elect_syst + govern_type2), # model 3
  update(m1_wbgee, .~. + elect_syst + govern_type2 + gini_disp), # model 4
  update(m1_wbgee, .~. + elect_syst + govern_type2 + ethn_fr), # model 5
  update(m1_wbgee, .~. + elect_syst + govern_type2 + gov_size_imf)) # model 6

### Table B12 appendix
modelsummary(wbre_models_no_authorit1, 
             stars = TRUE, 
             coef_map = coeff_labels_wbre,
             caption = "The effect of DPF on party corruption: Within-between random-effects models without authoritarian regimes (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             longtable = TRUE,
             output = "kableExtra") %>%
  kable_styling(font_size = 10, full_width = F, latex_options = "hold_position") %>%
  column_spec(column = 1, width = "4cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

###########################################
### Fixed-effects models lagged by BEEPS wave

### Fit lagged Year FE
models_lag_beeps <- list(
  m1_yfe_lag <- feols(capture_impact_1_norm ~ dpfv_dprc_infl_lag + private_agreg_lag + democracy_scaled_lag + log(gdppc_ppp_2017_lag) | year, data = impact_panel, cluster = ~country), ## Model 1
  
  m2_yfe_lag <- update(m1_yfe_lag, . ~ . + elect_syst_lag + govern_type2_lag + gov_size_imf_lag), # Model 2
  
  m3_cfe_lag <- feols(capture_impact_1_norm ~ dpfv_dprc_infl_lag + private_agreg_lag + democracy_scaled_lag + log(gdppc_ppp_2017_lag) | country, data = impact_panel, cluster = ~country), # Model 3
  
  m4_cfe_lag <- update(m3_cfe_lag, . ~ . + elect_syst_lag + govern_type2_lag + gov_size_imf_lag), # Model 4
  
  m5_tfe_lag <- feols(capture_impact_1_norm ~ dpfv_dprc_infl_lag + private_agreg_lag + democracy_scaled_lag + log(gdppc_ppp_2017_lag) | country + year, data = impact_panel, cluster = ~country), # Model 5
  
  m6_tfe_lag <- update(m5_tfe_lag, . ~ . + elect_syst_lag + govern_type2_lag + gov_size_imf_lag)) # Model 6

### Clean term labels
model_labels_lag_beeps <- c("dpfv_dprc_infl_lag" = "DPF level", 
                      "private_agreg_lag" = "Donation limits",
                      "democracy_scaled_lag" =  "Democracy",
                      "log(gdppc_ppp_2017_lag)" = "GDPpc(log)",
                      "elect_syst_lag" = "Electoral System",
                      "govern_type2_lag" = "Parliamentary", 
                      "ethn_fr_lag" = "Ethnic fraction.",
                      "gini_disp_lag" = "Inequality", 
                      "gov_size_imf_lag" = "Government size % GDP")

### Table B13 appendix
modelsummary(models_lag_beeps,
             stars = TRUE,
             coef_map = model_labels_lag_beeps,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: fixed-effects models lagged by BEEPS wave (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country.", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm") %>%
  add_header_above(c("", "Year fixed-effects" = 2, "Country fixed-effects" = 2, "Two-way fixed-effects" = 2))

############################################
### WBRE MODELS WITH PUBLIC FUNDING ADJUSTED FOR INFLATION, DISCOUNTED AND LAGGED BY BEEPS WAVES

#### Declaring panel data
panel_wbre_lag <- panel_data(impact, id = country, wave = year) 

### Fit WBRE lagged models

wbre_models_lag <- list(
  m1_wbgee <- wbgee(capture_impact_1_norm ~ dpfv_dprc_infl_lag + private_agreg_lag + democracy_scaled_lag + log(gdppc_ppp_2017_lag), data = panel_wbre_lag, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1"), # model 1
  update(m1_wbgee, .~. + elect_syst_lag), # model 2
  update(m1_wbgee, .~. + elect_syst_lag + govern_type2_lag), # model 3
  update(m1_wbgee, .~. + elect_syst_lag + govern_type2_lag + gini_disp_lag), # model 4
  update(m1_wbgee, .~. + elect_syst_lag + govern_type2_lag + ethn_fr_lag), # model 5
  update(m1_wbgee, .~. + elect_syst_lag + govern_type2_lag + gov_size_imf_lag))

### Clean term labels
coeff_labels_wbre_lag <- 
  c("dpfv_dprc_infl_lag" = "DPF level (within)",
    "private_agreg_lag" = "Donation limits (within)",
    "democracy_scaled_lag" =  "Democracy (within)", 
    "log(gdppc_ppp_2017_lag)" = "GDPpc(log) (within)",
    "elect_syst_lag" = "Electoral System (within)", 
    "govern_type2_lag" = "Parliamentary (within)",
    "ethn_fr_lag" = "Ethnic fractionalization (within)",
    "gini_disp_lag" = "Inequality (within)", 
    "gov_size_imf_lag" = "Government size % GDP (within)", 
    "imean(dpfv_dprc_infl_lag)" = "DPF level (between)", 
    "imean(private_agreg_lag)" = "Donation limits (between)",
    "imean(democracy_scaled_lag)" =  "Democracy (between)",
    "imean(log(gdppc_ppp_2017_lag))" = "GDPpc(log) (between)", 
    "imean(elect_syst_lag)" = "Electoral System (between)", 
    "imean(govern_type2_lag)" = "Parliamentary (between)", 
    "imean(ethn_fr_lag)" = "Ethnic fractionalization (between)", 
    "imean(gini_disp_lag)" = "Inequality (between)", 
    "imean(gov_size_imf_lag)" = "Government size % GDP (between)")

### Table B14 appendix

modelsummary(wbre_models_lag,
             stars = TRUE, 
             coef_map = coeff_labels_wbre_lag,
             coef_omit = "year|country|(Intercept)",
             caption = "The effect of DPF on party corruption: Within-between random-effects models lagged by BEEPS waves (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             longtable = TRUE,
             output = "kableExtra") %>%
  kable_styling(font_size = 10, full_width = F, latex_options = "hold_position") %>%
  column_spec(column = 1, width = "4cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

#############################################
### Fit fixed-effects models with three-years lagged and uniformly aggregated controls
### A new already processed dataset with 3-year lags
impact_lags_aggr <- read_dta("impact_lags_aggr.dta")
### Declare panel data
impact_panel_lags_aggr <- panel(impact_lags_aggr, ~country + year)

### Fit fixed-effects models with uniform lags
models_lags_aggr <- list(
  m1_yfe_lag_aggr <- feols(capture_impact_1_norm_lag ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | year, data = impact_panel_lags_aggr, cluster = ~ country), # Model 1
  
  m2_yfe_lag_aggr <- update(m1_yfe_lag_aggr, . ~ . + elect_syst + govern_type2 + gov_size_imf), # Model 2, 
  
  m3_cfe_lag_aggr <- feols(capture_impact_1_norm_lag ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country, data = impact_panel_lags_aggr, cluster = ~ country), # Model 3
  
  m4_cfe_lag_aggr <- update(m3_cfe_lag_aggr, . ~ . + elect_syst + govern_type2 + gov_size_imf), # Model 4
  
  m5_tfe_lag_aggr <- feols(capture_impact_1_norm_lag ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country + year, data = impact_panel_lags_aggr, cluster = ~ country), # Model 5
  
  m6_tfe_lag_aggr <- update(m5_tfe_lag_aggr, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 6

### Clean term labels
model_labels_lags_aggr <- 
  c("dpfv_dprc_infl" = "DPF level",
    "private_agreg" = "Donation limits", 
    "democracy_scaled" =  "Democracy", 
    "log(gdppc_ppp_2017)" = "GDPpc(log)",
    "elect_syst" = "Electoral System",
    "govern_type2" = "Parliamentary",
    "ethn_fr" = "Ethnic fraction.",
    "gini_disp" = "Inequality", 
    "gov_size_imf" = "Government size % GDP")

### Table B15 appendix
modelsummary(models_lags_aggr,
             stars = TRUE,
             coef_map = model_labels_lags_aggr,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: fixed-effects models with uniform lags (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country. Independent variable and controls are aggregated by uisng values from three years before each BEEPS wave", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "3cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm") %>%
  add_header_above(c("", "Year fixed-effects" = 2, "Country fixed-effects" = 2, "Two-way fixed-effects" = 2))

############################################
### Fit within-between random-effects models with three-years lagged and uniformly aggregated controls
### Declare panel data
panel_wbre_lags_aggr <- panel_data(impact_lags_aggr, id = country, wave = year) 

### Fit WBRE models with uniform lags
wbre_models_lags_aggr <- list(
  m1_wbgee_laggr <- wbgee(capture_impact_1_norm_lag ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017), data = panel_wbre_lags_aggr, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1"), # model 1
  update(m1_wbgee_laggr, .~. + elect_syst), # model 2
  update(m1_wbgee_laggr, .~. + elect_syst + govern_type2), # model 3
  update(m1_wbgee_laggr, .~. + elect_syst + govern_type2 + gini_disp), # model 4
  update(m1_wbgee_laggr, .~. + elect_syst + govern_type2 + gini_disp + ethn_fr), # model 5
  update(m1_wbgee_laggr, .~. + elect_syst + govern_type2 + gini_disp + gov_size_imf)) # model 6

### Clean term labels
coeff_labels_wbre_lags_aggr <- 
  c("dpfv_dprc_infl" = "DPF level (within)", 
    "private_agreg" = "Donation limits (within)",
    "democracy_scaled" =  "Democracy (within)", 
    "log(gdppc_ppp_2017)" = "GDPpc(log) (within)", 
    "elect_syst" = "Electoral System (within)", 
    "govern_type2" = "Parliamentary (within)", 
    "ethn_fr" = "Ethnic fractionalization (within)", 
    "gini_disp" = "Inequality (within)", 
    "gov_size_imf" = "Government size % GDP (within)",
    "imean(dpfv_dprc_infl)" = "DPF level (between)",
    "imean(private_agreg)" = "Donation limits (between)",
    "imean(democracy_scaled)" =  "Democracy (between)",
    "imean(log(gdppc_ppp_2017))" = "GDPpc(log) (between)",
    "imean(elect_syst)" = "Electoral System (between)",
    "imean(govern_type2)" = "Parliamentary (between)",
    "imean(ethn_fr)" = "Ethnic fractionalization (between)",
    "imean(gini_disp)" = "Inequality (between)",
    "imean(gov_size_imf)" = "Government size % GDP (between)")

### Table B16 appendix

modelsummary(wbre_models_lags_aggr, 
             stars = TRUE, 
             coef_map = coeff_labels_wbre_lags_aggr,
             caption = "The effect of DPF on party corruption: Within-between random-effects models with uniform lags (DPF inflation adjusted and discounted)",
             align = "lrrrrrr",
             longtable = TRUE,
             output = "kableExtra") %>%
  kable_styling(font_size = 10, full_width = F, latex_options = "hold_position") %>%
  column_spec(column = 1, width = "4cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  column_spec(column = 6, width = "1.7cm") %>%
  column_spec(column = 7, width = "1.7cm")

############################################
### Models with alternative operationalization of corruption, including, payments from central and local governmental officials
### Public funding inflation adjusted and discounted

### Fit TWFE models
models_twfe_altern <- list(
  m1_twfe_imp2 <- feols(capture_impact_2_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country + year, data = impact_panel, cluster = ~ country), # Model 1
  m2_twfe_imp2 <- update(m1_twfe_imp2, . ~ . + elect_syst + govern_type2 + gov_size_imf), # Model 2
  m3_twfe_imp3 <- feols(capture_impact_3_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) | country + year, data = impact_panel, cluster = ~ country),  # Model 3
  m4_twfe_imp3 <- update(m3_twfe_imp3, . ~ . + elect_syst + govern_type2 + gov_size_imf)) # Model 4

### Clean term labels
model_labels_alt <- 
  c("dpfv_nomin_infl" = "DPF level-inflation adjusted",
    "dpfv_dprc_infl" = "DPF level-inflation adjusted & discounted",
    "private_agreg" = "Donation limits",
    "democracy_scaled" =  "Democracy",
    "log(gdppc_ppp_2017)" = "GDPpc(log)",
    "elect_syst" = "Electoral System", 
    "govern_type2" = "Parliamentary",
    "ethn_fr" = "Ethnic fraction.",
    "gini_disp" = "Inequality",
    "gov_size_imf" = "Government size % GDP")

## Table B17 appendix
modelsummary(models_twfe_altern,
             stars = TRUE,
             coef_map = model_labels_alt,
             gof_omit = "Std.Error|R2 Within Adj",
             caption = "The effect of DPF on political corruption: two-way fixed-effects models with alternative operationalization of corruption",
             align = "lrrrr",
             output = "kableExtra") %>%
  kable_styling(font_size = 11, full_width = F, latex_options = "hold_position") %>%
  footnote(general = "Table entries represent unstandardised coefficients with robust standard errors clustered by country. Corruption 2 reflects the percentage of firms significantly and very significantly affected by the informal payments to political parties, legislators and central government officials. Corruption 3 reflects the percentage of firms significantly and very significantly affected by the informal payments to political parties, legislators, central and local government officials", footnote_as_chunk = TRUE, threeparttable = T) %>%
  column_spec(column = 1, width = "6cm") %>%
  column_spec(column = 2, width = "1.8cm") %>% 
  column_spec(column = 3, width = "1.8cm") %>% 
  column_spec(column = 4, width = "1.8cm") %>% 
  column_spec(column = 5, width = "1.8cm") %>%
  add_header_above(c("", "Corruption 2" = 2, "Corruption 3" = 2))

############################################
## Within-Between Random Effects Models (WBRE) with alternative operationalization of corruption 
### Public funding inflation adjusted only

# Fit WBRE 
wbre_models_altern <- list(
  m1_wbgee_imp2 <- wbgee(capture_impact_2_norm ~ dpfv_nomin_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017), data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1"), # model 1
  m2_wbgee_imp2 <- update(m1_wbgee_imp2, .~. + elect_syst + govern_type2 + gov_size_imf), # model 2
  m3_wbgee_imp3 <- wbgee(capture_impact_3_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017), data = panel_wbre, model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1"), # model 3
  m4_wbgee_imp3 <- update(m3_wbgee_imp3, .~. + elect_syst + govern_type2 + gov_size_imf)) # model 4

# Clean term labels 
coeff_labels_wbgee_altern <- 
  c("dpfv_nomin_infl" = "DPF level-inflation adjusted (within)",
    "dpfv_dprc_infl" = "DPF level-inflation adjusted & discounted (within)", 
    "private_agreg" = "Donation limits (within)", 
    "democracy_scaled" =  "Democracy (within)", 
    "log(gdppc_ppp_2017)" = "GDPpc(log) (within)", 
    "elect_syst" = "Electoral System (within)", 
    "govern_type2" = "Parliamentary (within)", 
    "gov_size_imf" = "Government size % GDP (within)", 
    "imean(dpfv_nomin_infl)" = "DPF level-inflation adjusted (between)",
    "imean(dpfv_dprc_infl)" = "DPF level-inflation adjusted & discounted (between)", 
    "imean(private_agreg)" = "Donation limits (between)",
    "imean(democracy_scaled)" =  "Democracy (between)", 
    "imean(log(gdppc_ppp_2017))" = "GDPpc(log) (between)", 
    "imean(elect_syst)" = "Electoral System (between)", 
    "imean(govern_type2)" = "Parliamentary (between)", 
    "imean(gov_size_imf)" = "Government size % GDP (between)")

### Table B18 appendix
modelsummary(wbre_models_altern, 
             stars = TRUE, 
             coef_map = coeff_labels_wbgee_altern,
             caption = "The effect of DPF on political corruption: within-between random-effects models with alternative operationalization of corruption",
             align = "lrrrr",
             longtable = TRUE,
             notes = list("Note: Corruption operationalization as for table B17."),
             output = "kableExtra") %>%
  kable_styling(font_size = 10, full_width = F, latex_options = "hold_position") %>%
  column_spec(column = 1, width = "6cm") %>%
  column_spec(column = 2, width = "1.7cm") %>% 
  column_spec(column = 3, width = "1.7cm") %>% 
  column_spec(column = 4, width = "1.7cm") %>% 
  column_spec(column = 5, width = "1.7cm") %>%
  add_header_above(c("", "Corruption 2" = 2, "Corruption 3" = 2))

############################################
### Fit country-based Jackknife estimations 

### ROBUSTNESS/SENSITIVITY TESTS for year FE models; country-based Jackknife

models_year_fe <- list() #list of models 
### For loop to iterate across all countries by dropping out a country at a time 
for(country_iterator in unique(impact_panel$country)) { #for loop leaving out all obs from a given country
  models_year_fe[[paste0("Drop ", country_iterator)]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | year, data = impact_panel %>% filter(country != country_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for plotting
models_year_fe_df <- models_year_fe %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep only public funding estimates
  mutate(period = "A: Year fixed-effects models", ### model type
         main_estim = -0.20) ### main analysis public funding estimate

#########################################
### ROBUSTNESS/SENSITIVITY TESTS for country FE models, country-based Jackknife

models_country_fe <- list() #list of models 
### Loop to iterate across all countries by dropping out a country at a time 
for(country_iterator in unique(impact_panel$country)) { #for loop leaving out all obs from a given country
  models_country_fe[[paste0("Drop ", country_iterator)]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | country, data = impact_panel %>% filter(country != country_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for plotting
models_country_fe_df <- models_country_fe %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep only public funding estimates
  mutate(period = "B: Country fixed-effects models", ### model type
         main_estim = -0.107) ### main analysis public funding estimate

#########################################
### ROBUSTNESS/SENSITIVITY TESTS for TWFE models, country-based Jackknife

models_twoway_fe <- list() #list of models 
### Loop to iterate across all countries by dropping out a country at a time
for(country_iterator in unique(impact_panel$country)) { #for loop leaving out all obs from a given country
  models_twoway_fe[[paste0("Drop ", country_iterator)]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | country + year, data = impact_panel %>% filter(country != country_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for plotting
models_twoway_fe_df <- models_twoway_fe %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep only public funding estimates
  mutate(period = "C: Two-way fixed-effects models", ### model type
         main_estim = -0.145) ### main analysis public funding estimate

###################

### Joining the data frames from the country-based Jackknife estimations for plotting

models_all_fe_df <- 
  bind_rows(models_year_fe_df,
            models_country_fe_df, 
            models_twoway_fe_df) %>%
  mutate(conf_low_90 = estimate - 1.645*std.error, # create 90% CI
         conf_high_90 = estimate + 1.645*std.error) # create 90% CI

### Plot the coefficient estimates for all models with a discarded country at a time, Figure C1 in Appendix - all fixed-effects models

ggplot(models_all_fe_df, aes(x = fct_reorder(id, estimate), y= estimate)) +
  geom_point(color = "grey40", size = 2, position = position_dodge(width = 1)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0), color = "grey40",  linewidth = 0.4, position = position_dodge(width = 1)) +
  geom_errorbar(aes(ymin = conf_low_90, ymax = conf_high_90, width = 0), color = "grey40",  linewidth = 1, position = position_dodge(width = 1)) +
  geom_hline(aes(yintercept = main_estim), linetype = 2, color = "navy") +
  geom_hline(yintercept = 0, linetype = 5, color = "red4") +
  facet_wrap(~ period, scales = "free_y", nrow = 3) +
  theme_minimal(base_family = "serif", base_size = 11) +
  theme(panel.border = element_rect(color = "gray30", fill = NA),
        panel.grid = element_blank(),
        strip.text = element_text(size = 11, color = "gray10"),
        plot.caption.position = "plot",
        plot.caption = element_text(size = 11, hjust = 0),
        axis.text.y = element_text(size = 11, hjust = 0),
        axis.text.x = element_text(size = 11, angle = 90),
        axis.title.y = element_blank()) +
  labs(x = "", y = "Public funding estimates with 90% and 95% Confidence Intervals",
       caption = str_wrap("Note: The blue line represents the public funding estimate from the main analysis and Appendix B (Table B1, Model 6 - Year fixed-effects, Table B2, Model 6 - Country fixed-effects, Table 1, Model 4 - Two-way fixed-effects). The red line indicates statistical significance for each model when the respective country is discarded. Figure entries represent unstandardised coefficients with robust standard errors clustered by country.", width = 120))

###########################################
### Fit observation-based Jackknife estimations

### ROBUSTNESS/SENSITIVITY TESTS for year FE models, observation-based Jackknife

jackknife_yfe_simple <- list() #list of models 

### Loop to iterate across all observations by dropping out an observation at a time
for(observation_iterator in unique(impact_panel$unique_id)) { #for loop leaving out one observation at a time
  jackknife_yfe_simple[[observation_iterator]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | year, data = impact_panel %>% filter(unique_id != observation_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for plotting
jackknife_yfe_simple_df <- jackknife_yfe_simple %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep public funding
  mutate(period = "A: Year fixed-effects models", ### Model type 
         main_estim = -0.20, ### main analysis public funding estimate
         conf_low_90 = estimate - 1.645*std.error, ### estimate 90% CI
         conf_high_90 = estimate + 1.645*std.error, ### estimate 90% CI
         id = as.numeric(id))

#########################################
### ROBUSTNESS/SENSITIVITY TESTS for country FE models, observation-based Jackknife

jackknife_cfe_simple <- list() #list of models 

### Loop to iterate across all observations by dropping out an observation at a time
for(observation_iterator in unique(impact_panel$unique_id)) { #for loop leaving out one observation at a time
  jackknife_cfe_simple[[observation_iterator]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | country, data = impact_panel %>% filter(unique_id != observation_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for plotting
jackknife_cfe_simple_df <- jackknife_cfe_simple %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep public funding
  mutate(period = "B: Country fixed-effects models", ### Model type
         main_estim = -0.107, ### main analysis public funding
         conf_low_90 = estimate - 1.645*std.error, ### estimate 90% CI
         conf_high_90 = estimate + 1.645*std.error, ### estimate 90% CI
         id = as.numeric(id))

##############################################
### ROBUSTNESS/SENSITIVITY TESTS for country TWFE models, observation-based Jackknife

jackknife_tfe_simple <- list() #list of models 

### Loop to iterate across all countries by dropping out an observation at a time
for(observation_iterator in unique(impact_panel$unique_id)) { #for loop leaving out one observation at a time
  jackknife_tfe_simple[[observation_iterator]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | country + year, data = impact_panel %>% filter(unique_id != observation_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for plotting
jackknife_tfe_simple_df <- jackknife_tfe_simple %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep public funding
  mutate(period = "C: Two-way fixed-effects models", ### Model type
         main_estim = -0.145, ### main analysis public funding
         conf_low_90 = estimate - 1.645*std.error, ### estimate 90% CI
         conf_high_90 = estimate + 1.645*std.error, ### estimate 90% CI
         id = as.numeric(id))

###########################
### Joining the data frames from the observation-based Jackknife estimations for plotting
jackknife_all_fe_df <- bind_rows(jackknife_yfe_simple_df,
                                 jackknife_cfe_simple_df,
                                 jackknife_tfe_simple_df) 
### Plot the coefficient estimates for all 130 models with a discarded observation at a time, Figure C2 in Appendix - all fixed-effects models

ggplot(jackknife_all_fe_df, aes(x = id, y = estimate)) +
  geom_ribbon(aes(ymin = conf.low,  ymax = conf.high), fill = "grey90", alpha = 0.5) +
  geom_ribbon(aes(ymin = conf_low_90,  ymax = conf_high_90), fill = "grey80", alpha = 0.5) +
  scale_x_continuous(breaks = seq(0, 130, 20)) +
  geom_line(size = 1, color = "grey20") +
  geom_hline(aes(yintercept = main_estim), linetype = 2, color = "navy") +
  geom_hline(aes(yintercept = 0), color = "red4", linetype = 5) +
  facet_wrap(~ period, scales = "free_y", ncol = 1) +
  theme_minimal(base_family = "serif") +
  theme(panel.border = element_rect(color = "gray30", fill = NA),
        panel.grid = element_blank(),
        strip.text = element_text(size = 12, color = "gray10"),
        plot.caption.position = "plot",
        plot.caption = element_text(size = 11, hjust = 0),) +
  labs(x = "Number of estimations", y = "Public funding estimates with 90% and 95% confidence intervals",
       caption = str_wrap("Note: The dark-grey line shows the variation in public funding coeffcients when each subsequent observation is dropped from the estimation. The dark-grey and light-grey shaded areas represent the confidence intervals at 90% and 95% levels respectively. The blue line indicates the estimates from the main analysis, while the the red line shows statistical significance for each model when every subsequent observation is discarded. The plot is based on Table B1, Model 6 - Year fixed-effects, Table B2, Model 6 - Country fixed-effects, Table 1, Model 4 - Two-way fixed-effects", width = 120))

#######################################
### ROBUSTNESS/SENSITIVITY TESTS FOR WITHIN-BETWEEN MODELS, country-based Jackknife

models_wbre <- list() # list of models
### Loop to iterate across all countries by dropping out a country at a time
for(country_iterator in unique(panel_wbre$country)) {#for loop leaving out all observations from a given country
  models_wbre[[ paste0("Drop ",country_iterator) ]] <- wbgee(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf, data = panel_wbre %>% filter(country != country_iterator), model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1")
}

### Extracting the public funding estimates into data frame for plotting
models_wbre_df <- lapply(models_wbre, '[[', 'coefs') %>% 
  bind_rows(.id = "id") %>% 
  mutate(main_estim_b = -0.316, # dpf estimate main analysis between effect
         main_estim_w = -0.107, # dpf estimate main analysis within effect
         term = as.character(rownames(.))) %>% 
  select(id, term, everything()) %>%
  filter(term != '(Intercept)', str_detect(term, "dpfv_dprc_infl")) %>% 
  mutate(term = str_remove_all(term, "\\.{3}\\d{1,}|`"),
         estim_type = ifelse(str_detect(term, "imean"), "Between", "Within"),
         term = ifelse(str_detect(term, "imean"), "DPF level (between)", "DPF level (within)"),
         conf_low = Est. - 1.96*S.E., ### estimate 95% CI
         conf_high = Est. + 1.96*S.E.,### estimate 95% CI
         conf_low_90 = Est. - 1.645*S.E., ### estimate 90% CI
         conf_high_90 = Est. + 1.645*S.E.) ### estimate 90% CI

### Plot the coefficient estimates for all models with a discarded country at a time, Figure C3, Panel A in Appendix - WBRE

p1_wbre <- models_wbre_df %>%
  filter(estim_type == "Between") %>%
  ggplot(aes(x = fct_reorder(id, Est.), y= Est.)) +
  geom_point(colour = "grey40", size = 2) +
  geom_errorbar(aes(ymin = conf_low, ymax = conf_high, width = 0), color = "grey40", size = 0.4) +
  geom_errorbar(aes(ymin = conf_low_90, ymax = conf_high_90, width = 0, color = estim_type), size = 0.8, position = position_dodge(width = 0.6)) +
  scale_color_grey(start = 0.4, end = 0.4) +
  scale_shape_discrete() + 
  geom_hline(aes(yintercept = main_estim_b), linetype = 2, color = "navy") +
  geom_hline(yintercept = 0, linetype = 5, color = "red4") +
  theme_minimal(base_family = "serif", base_size = 11) +
  theme(panel.border = element_rect(color = "gray30", fill = NA),
        panel.grid = element_blank(),
        strip.text = element_text(size = 11, color = "gray10"),
        strip.background = element_rect(fill = "white"),
        legend.position = "none",
        plot.caption.position = "plot",
        plot.caption = element_text(size = 11, hjust = 0),
        axis.text.y = element_text(size = 11),
        axis.text.x = element_text(size = 11, vjust = 0.5, angle = 90)) +
  labs(title = "A: Country-based Jakknife variance estimation (Between effect)", x = "" , y = "Public funding estimates with \n90% and 95% Confidence Intervals", 
       caption = str_wrap("Note: The blue line represents the public funding estimate from the main analysis (Table 2, Model 4). The red line indicates statistical significance for each model when the respective country is discarded. Figure entries represent unstandardised coefficients with robust standard errors clustered by country.", width = 120))

###################################
### ROBUSTNESS/SENSITIVITY TESTS FOR WITHIN-BETWEEN MODELS, observation-based Jackknife

models_wbre_simple <- list() # list of models
### Loop to iterate across all observations by dropping out an observation at a time
for(observation_iterator in unique(panel_wbre$unique_id)) {#for loop leaving out an observation at a time
  models_wbre_simple[[observation_iterator]] <- wbgee(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf, data = panel_wbre %>% filter(unique_id != observation_iterator), model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1")
}

### Extracting the public funding estimates into data frame for plotting
models_wbre_simple_df <- lapply(models_wbre_simple, '[[', 'coefs') %>% 
  bind_rows(.id = "id") %>% 
  mutate(main_estim_b = -0.316, # dpf estimate main analysis between effect
         main_estim_w = -0.107, # dpf estimate main analysis within effect, not shown in figure since it is equivalent of country FE
         term = as.character(rownames(.))) %>% 
  select(id, term, everything()) %>%
  filter(term != '(Intercept)', str_detect(term, "dpfv_dprc_infl")) %>% 
  mutate(id = as.numeric(id),
         term = str_remove_all(term, "\\.{3}\\d{1,}|`"),
         estim_type = ifelse(str_detect(term, "imean"), "Between", "Within"),
         term = ifelse(str_detect(term, "imean"), "DPF level (between)", "DPF level (within)"),
         conf_low = Est. - 1.96*S.E., ### estimate 95% CI
         conf_high = Est. + 1.96*S.E.,### estimate 95% CI
         conf_low_90 = Est. - 1.645*S.E.,### estimate 90% CI
         conf_high_90 = Est. + 1.645*S.E.) %>% ### estimate 90% CI
  filter(estim_type == "Between")

### Plot the coefficient estimates for all models with a discarded observation at a time, Figure C3, panel B in Appendix - WBRE

p2_wbre <- ggplot(models_wbre_simple_df, aes(x = id, y = Est.)) +
  geom_ribbon(aes(ymin = conf_low,  ymax = conf_high), fill = "grey90", alpha = 0.5) +
  geom_ribbon(aes(ymin = conf_low_90,  ymax = conf_high_90), fill = "grey80", alpha = 0.5) +
  scale_x_continuous(breaks = seq(0, 130, 20)) +
  geom_line(size = 1, color = "grey20") +
  geom_hline(aes(yintercept = main_estim_b), linetype = 2, color = "navy") +
  geom_hline(aes(yintercept = 0), color = "red4", linetype = 5) +
  theme_minimal(base_family = "serif") +
  theme(panel.border = element_rect(color = "gray30", fill = NA),
        panel.grid = element_blank(),
        strip.text = element_text(size = 12, color = "gray10"),
        plot.caption.position = "plot",
        plot.caption = element_text(size = 11, hjust = 0)) +
  labs(title = "B: Observation-based Jakknife variance estimation (Between effect)", x = "Number of estimations", y = "Public funding estimates with \n90% and 95% Confidence Intervals",
       caption = str_wrap("Note: The dark-grey line shows the variation in public funding coeffcients when each subsequent observation is dropped from the estimation. The dark-grey and light-grey shaded areas represent the confidence intervals at 90% and 95% levels respectively. The blue line represents the public funding estimate from the main analysis (Table 2, Model 4) The red line indicates statistical significance for each model when every subsequent observation is discarded.", width = 120))

### Combine WBRE Jackknife plots 
plot_grid(p1_wbre, p2_wbre, nrow = 2)

##########################################
### ROBUSTNESS/ SENSITIVITY TESTS for BEEPS round-based Jackknife estimations

### Year FE
year_fe_drop_beeps_rounds <- list() #list of models 

### Loop to iterate across all BEEPS rounds by dropping out a BEEPS round at a time 
for(year_iterator in unique(impact_panel$year)) { #for loop leaving out all obs from a given country
  year_fe_drop_beeps_rounds[[paste0("Drop ", year_iterator)]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | year, data = impact_panel %>% filter(year != year_iterator), cluster = ~country))
}
### Extracting the public funding estimates into data frame for plotting
year_fe_drop_beeps_rounds_df <- year_fe_drop_beeps_rounds %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep only public funding estimates
  mutate(estim_type = "B: Year fixed-effects", ### model type
         main_estim = -0.20) ### main analysis public funding estimate

######################################
### Country FE
country_fe_drop_beeps_rounds <- list() #list of models 

### Loop to iterate across all BEEPS rounds by dropping out a BEEPS round at a time 
for(year_iterator in unique(impact_panel$year)) { #for loop leaving out all obs from a given country
  country_fe_drop_beeps_rounds[[paste0("Drop ", year_iterator)]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | country, data = impact_panel %>% filter(year != year_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for plotting
country_fe_drop_beeps_rounds_df <- country_fe_drop_beeps_rounds %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep only public funding estimates
  mutate(estim_type = "A: Country fixed-effects", ### model type
         main_estim = -0.107) ### main analysis public funding estimate

##########################################
### TWFE

twfe_fe_drop_beeps_rounds <- list() #list of models 
### Loop to iterate across all countries by dropping out a country at a time full data.
for(year_iterator in unique(impact_panel$year)) { #for loop leaving out all obs from a given country
  twfe_fe_drop_beeps_rounds[[paste0("Drop ", year_iterator)]] <- tidy(conf.int = TRUE, feols(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf | country + year, data = impact_panel %>% filter(year != year_iterator), cluster = ~country))
}

### Extracting the public funding estimates into data frame for further analysis and plotting
twfe_fe_drop_beeps_rounds_df <- twfe_fe_drop_beeps_rounds %>%
  bind_rows(.id = "id") %>%
  filter(term == "dpfv_dprc_infl") %>% ### keep only public funding estimates
  mutate(estim_type = "C: Two-way fixed-effects", ### model type
         main_estim = -0.145) ### main analysis public funding estimate

####################################
### ROBUSTNESS/ SENSITIVITY TESTS for WBRE models, BEEPS round-based Jackknife 
models_wbre_drop_beeps <- list() # list of models

for(year_iterator in unique(panel_wbre$year)) {#for loop leaving out all observations from a given country
  models_wbre_drop_beeps[[ paste0("Drop ", year_iterator) ]] <- wbgee(capture_impact_1_norm ~ dpfv_dprc_infl + private_agreg + democracy_scaled + log(gdppc_ppp_2017) + elect_syst + govern_type2 + gov_size_imf, data = panel_wbre %>% filter(year != year_iterator), model = "w-b", use.wave = F, scale = FALSE, balance.correction = TRUE, cor.str = "ar1")
}

### Data manipulation to extract public funding estimates 
models_wbre_drop_beeps_df <- lapply(models_wbre_drop_beeps, '[[', 'coefs') %>% 
  bind_rows(.id = "id") %>% 
  mutate(term = as.character(rownames(.))) %>% 
  select(id, term, everything()) %>%
  filter(term != '(Intercept)', str_detect(term, "dpfv_dprc_infl")) %>% 
  mutate(term = str_remove_all(term, "\\.{3}\\d{1,}|`"),
         estim_type = ifelse(str_detect(term, "imean"), "Between", "Within"),
         term = ifelse(str_detect(term, "imean"), "DPF level (between)", "DPF level (within)"),
         conf.low = Est. - 1.96*S.E., ### estimate confidence intervals
         conf.high = Est. + 1.96*S.E., ### estimate confidence intervals
         main_estim = -0.316) %>% ### dpf estimate from main analysis estimate
  filter(estim_type == "Between") %>% ### Effect type
  mutate(estim_type = "D: Within-bewteen random-effects") %>% ### Model type
  rename("estimate" = "Est.", "std.error" = "S.E.", "statistic" = "z val.", "p.value" = "p") ### rename variables
### clean row names
row.names(models_wbre_drop_beeps_df) <- c(1,2,3,4,5) 

###############################
### Joining the data frames for all models BEEPS round-based Jackknife estimations - year FE, country FE, TWFE, WBRE

models_all_drop_beeps <- 
  bind_rows(year_fe_drop_beeps_rounds_df, 
            country_fe_drop_beeps_rounds_df, 
            twfe_fe_drop_beeps_rounds_df,
            models_wbre_drop_beeps_df) %>%
  mutate(conf_low_90 = estimate - 1.645*std.error, ### create 90% conf int
         conf_high_90 = estimate + 1.645*std.error) ### create 90% conf int

### Plotting the coefficient estimates for all models with a discarded BEEPS round at a time, Figure C4 in Appendix 

ggplot(models_all_drop_beeps, aes(x = fct_rev(id), y= estimate)) +
  geom_point(color = "grey40", size = 2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = 0), color = "grey40",  size = 0.4) +
  geom_errorbar(aes(ymin = conf_low_90, ymax = conf_high_90, width = 0), color = "grey40",  size = 1) +
  geom_hline(aes(yintercept = main_estim), linetype = 2, color = "navy") +
  geom_hline(yintercept = 0, linetype = 5, color = "red4") +
  coord_flip() +
  facet_wrap(~ estim_type, ncol = 2) +
  theme_light(base_family = "serif") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, color = "gray10", face = "bold"),
        plot.caption.position = "plot",
        plot.caption = element_text(size = 11, hjust = 0),
        axis.text = element_text(size = 12, hjust = 0.5, color = "gray10"),
        axis.title.x = element_text(size = 12, hjust = 0.5),
        axis.title.y = element_blank()) +
  labs(x = "", y = "Public funding estimates with 90% and 95% Confidence Intervals",
       caption = str_wrap("Note: Jackknife variance estimations based on BEEPS rounds. The blue line represents the public funding estimate from the Appendix B and main analysis respectively (Table B1, Model 6 - Year fixed-effects; Table B2, Model 6 - Country fixed-effects; Table 1, Model 4 (main text) - Two-way fixed-effects; Table 2, Model 4 (main text) - Within-between random-effects). The red line indicates wether the estimates are statistically significant for each model when the respective BEEPS round is removed from the estimation. Figure entries represent unstandardised coefficients with robust standard errors clustered by country.", width = 120))

